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F^ Abstract 

' We present a methodology to change the state of the Weather Research Forecasting (WRF) model coupled with 

b^ the fire spread code SPIRE, based on Rothermel's formula and the level set method, and with a fuel moisture model. 

• The fire perimeter in the model changes in response to data while the model is running. However, the atmosphere 

Q state takes time to develop in response to the forcing by the heat flux from the fire. Therefore, an artificial fire history 

' ^ is created from an earlier fire perimeter to the new perimeter, and replayed with the proper heat fluxes to allow the 

^^ atmosphere state to adjust. The method is an extension of an earlier method to start the coupled fire model from a 

F^ developed fire perimeter rather than an ignition point. The level set method can be also used to identify parameters 

of the simulation, such as the fire spread rate. The coupled model is available from |openwf m. org| and it extends the 

WRF-Fire code in WRF release. 
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O 1. Introduction 

(N 

^^ This article reports on recent developments in building a Dynamic Data Driven Application System (DDDAS) for 

^ wildland fire simulations |[T]|2l|3l. A DDDAS is based on the ability to incorporate data into an executing simulation 

H. See Fig.[T]for the overaU scheme of the DDDAS. 

The paper is organized as follows. In Sec. [2j we we review some existing approaches to data assimilation in 

C^ simulations of wildland fires. In Sec. [3]and|4] we briefly formulate the model and the principal idea of creating 

and replaying artificial fire history from point ignition to a given perimeter, for reference. Sec. |5] considers several 

methods for the construction of a level set function, needed for the replay, two from our previous work and a new 

method, which is an extension of the reinitialization equation approach known in level set methods. In Sec. [6j we 

present a new method how the level set functions constructed for two perimeters can be used to create and replay 
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METEO INPUT DATA 



Large scale weather data from 
NOAA Rapid Update Cycle runs: 

• 12km-resolution initial conditions 

• 12km-resolution boundary 
conditions 

Static data: 

• High-resolution topography 

• Land Use and Soil Data 
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WRF SFIRE 
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• High resolution fuel data: 

• 30m-resolution fuel description 

• 30m-resolution elevation data 
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DYNAMIC FIRE DATA 



WRF framework (atmosphere) 
ARW atmospheric core 
WPS preprocessing system 



4'KFIRE-AFFECTED| 
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|f ire-g enerated * 
iheat and moisture 



Fire Model: 
Rothermel fire spread model 





Fire front trackina based on 
the level set metnod 



• Satellite & aircraft fire detection 

• NDWI-derived fuel moisture, NDVI- 
corrected fuel data 

• MiSR plume heights and fire 
intensity 

• Incident fire perimeters 

• Ground sensors 



r-C 



DATA ASSIMILATION 



• Position correction 

• Ensemble Kalman filters 

• Wavelet and spectral filters 

• Perimeter adjustment 



METEO OUTPUT 



High-resolution forecast: 

• wind speed and direction 

• air temperature 

• air humidity 

• precipitation 

• cloudiness etc... 



High-resolution fire forecast: 
fire area 
fire risk 
fire heat flux 
fire intensity 
fire rate of spread 
Diume heiaht 



Figure 1: Scheme of wildland fire DDDAS. 



an artificial fire history between the two perimeters, and a new method that uses the two level set functions for an 
automatic adjustment of the fire spread rate between the two perimeters. Sec. |7] describes a new moisture model 
coupled with the fire and atmosphere model, and the possibilities for the assimilation of moisture data. Finally, Sec. [8] 
is the conclusion. 



2. Data assimilation for wildland fires 

One way to incorporate data into an executing simulation is by sequential statistical estimation, which takes all 
available data to date into account, and is known in geosciences as data assimilation. Data assimilation is a standard 
technique in numerical weather prediction, and the ability to assimilate large amounts of real-time data is behind 
much of the recent improvement of weather forecast skill O. However, data assimilation for wildland fires poses 
unique challenges, and classical data assimilation methods simply will not work (61121 [3. One of the reasons is that in 
many other physical systems where standard methods work well, such as pollution transport or atmospheric dynamics, 
unwanted perturbations tend to dissipate over time; but, in a fire model, once a perturbation ignites an unwanted fire, 
the fire will keep growing, and after few assimilation cycles, everything burns. Another reason is that a fire as a 
coherent structure, needs to be moved, started, or extinguished in response to the data, which requires positional, 
Lagrangean correction; additive corrections of the values of the physical fields are not very useful. 

Data assimilation methods by sequential Monte-Carlo methods (SMC), also known as particle filters, were 
developed in the literature for cell-based fire models ||9l[T0l. They can handle non-Gaussian distributions, but they 
are computationally very expensive, because they require very large ensembles to cover a region of the state space 
by random perturbations. A suitable perturbation algorithm is the key to a successful application. The perturbation 
methods used in wildland fire modeling range from random modifications of the burn area 1 9 1 to genetic algorithms, 
which evolve the shape of the fire by simulated evolution, where the states with fire regions closer the the data are more 
likely to survive [11 J. While SMC methods with tens of thousands of particles may be feasible for 2D cell models, 
with relatively small state vectors, they are definitely out of question for a coupled atmosphere-fire model. Methods 
based on the optimal statistical interpolation and the Kalman filter (KF), such as the ensemble Kalman filter (EnKF), 
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Figure 2: Visualization in Google Earth client, 2007 Witch Fire. False color shows fire heat flux, superimposed on the 
Earth surface. Patches of slower fuels keep burning behind the fireline. Reproduced from ifTSl . 



assume that the state distribution is at least approximately Gaussian and they modify the state in response to data 
(3 p. 180] rather than rely on hitting the right answer with random perturbations. Thus, KF-based mehods require 
much smaller ensembles that SMC methods, but still in the range of 20-100 members and easily many hundreds 
lfT2ll . However, because of the fine resolution of the atmospheric model needed over large areas, and the associated 
need for small time steps, the simulations are computationally very demanding, and such ensembles are still out of 
question. FFT-based data assimilation methods, which reduce data assimilation to efficient operations with diagonal 
matrices |[T3ll and can drastically reduce the required ensemble size, from hundreds to often just 5 or 10 members. 
However, using the Fourier basis is tantamount to the assumption that the state co variance does not vary spatially 1 14|. 
Wavelet estimation can combine the eff'ectiveness of spectral methods with an automatic treatment of spatial locality 
ifTSll . Wavelet diagonal approximations of the covariance matrix |[T6ll are of particular interest, as they allow efficient 
evaluation of the EnKF formulas 1 17|. 

Position correction methods, such as morphing |6|, can overcome the limitations of changing the state of the 
simulation by additive corrections only. These method extend the state by a new variable containing a deformation 
field, similarly as in optical flow methods 1 19] and extraction of the wind field from a sequence of radar images 
1201 . For other related position correction methods, see, e.g., 1211 . Our morphing technique is distinguished by 
replacing linear combinations of member states, which are at the heart of, e.g., the EnKF, by intermediate states, 
which interpolate both the magnitude and the position of coherent features, such as fires. Time series of station 
observations could be handled by considering composite states over several time steps. However, while morphing 
works successfully for fire models (HHl, it changes the delicate physical balance of the atmospheric equations and 
limits the possibility of the treatment of the model as a black box. Even a simple linear transformation to move 
and reshape the vortex in hurricane forecasting needs rebalancing of the atmospheric variables from conservation 
equations |[22l p. 11]. 

Therefore, an important problem in data assimilation for a coupled atmosphere-fire model is how to adjust the 
atmosphere state when the state of the fire model changes in response to data. The heat output of the fire is concentrated 
in a narrow area with active combustion, therefore the fire forcing on the atmosphere is highly localized. If the fire 
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is just shifted, a position correction alone can be successful! to some extent (61 [H because the relationship between 
the changes in the atmosphere and in the fire is captured in the covariance of their deformation fields. However, in 
general, the covariance does not contain sufficient information and a spin-up is required to develop proper circulation 
patterns for the changed fire forcing. 

3. The coupled atmosphere - fire model 

Over time, the wildland fire DDDAS has evolved from a simple convection-reaction-diff'usion equation exploratory 
model to test data assimilation methodologies I.23J and the CAWFE model (241 [25], which couples the Clark-Hall 
atmospheric model with fire spread implemented by tracers (Lagrangean particles), to the currently used Weather 
Research Forecasting (WRF) mesoscale atmospheric code 1 26] coupled with a spread model implemented by the level 
set method L8,,27J. The Clark-Hall model has many favorable properties, such as the ability to handle refinement, but 
WRF is a supported community model, it can execute in parallel, and has built-in export and import of state, which is 
essential for data assimilation. Also, WRF supports data formats standard in geosciences. The implementation by the 
level set method was chosen because the level set function can be manipulated much more easily than tracers. The 
coupled code is available from the Open Wildland Fire Modeling Environment (OpenWFM) (TSl at openwf m . org} 
which contains also diagnostic and data processing utilities, including visualization in Google Earth (Fig. |2]), which 
we first proposed in |2|. A subset of the SFIRE code was released with WRF as WRF-Fire. The model is capable of 
running on a cluster faster than real time with atmospheric resolution in tens of m, needed to resolve the atmosphere- 
fire interaction, for a fire of size over 10 km [28 1. See (8l[23 for futher details and references. 

The state variables of the fire model are the level set function, O, the time of ignition Ti, and the fuel fraction 
remaining F, given by their values on the nodes of the fire model mesh. At a given simulation time t, the fire area is 
represented by the level set function O as the set of all points (x, y) where O (t, x, y) < 0. Since the level set function 
is interpolated linearly between nodes, this allows a submesh representation of the fire area. In every time step of the 
simulation, the level set function is advanced by one step of a Runge-Kutta scheme for the level set equation 

JO 

-- = -7?|VO|, (1) 

at 

where R = R(t, x,y) is the fire rate of spread and H is the Euclidean norm. The ignition time Ti = T^ (x, j) is then 
computed for all newly ignited nodes, and it satisfies the consistency condition 

O a, x,y)<0<^^Ti {x,y) < t, (2) 

where both inequalities express the condition that the location (x, y) is burning at the time t. 

The fire rate of spread R is given by the Rothermel's formula |29 1 as a function of the wind speed (at a height 
dependent on the fuel) and the slope in the direction normal to the fireline. From the level-set representation of the 
fireline at the time ^ as O (t, x,y) = 0, it follows by an easy calculus that the normal direction is VO/ |VO|, where H is 
the Euclidean norm. Thus, 

where u is the wind field. 

Once the fuel starts burning, the remaining mass fraction F = F(t, x, y) is approximated by exponential decay, 

^ I 1, t<Ti{x,y), 

where Tf is the fuel burn time, i.e., the number of seconds for the fuel to burn down to 1/e ^ 0.3689 of the starting 
fuel fraction F = I. The heat fluxes from the fire to the atmosphere are taken proportional to the fuel burning rate, 
dF (t, x,y) /dt. The proportionality constants are fuel coefficients. The heat fluxes from the fire are inserted into the 
atmospheric model as forcing terms in diff'erential equations of the atmospheric model in a layer above the surface, 
with exponential decay with altitude. This scheme is required because atmospheric models with explicit timestepping, 
such as WRF, do not support flux boundary conditions. The sensible heat flux is added to the time derivative of the 
temperature, while the latent heat flux is added to the derivative of water vapor concentration. 
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Figure 3: Creating artificial time history by proportions. The ignition times at the ignition point A and the fire 
perimeter T are interpolated linearly along the segment between A and a point 5 on F to a mesh point X. 

4. Replaying artificial fire history 

The SFIRE code as presented in IH [27l starts from one or more ignition points. The release of the heat from 
the fire then gradually establishes atmospheric circulation patterns and the fire evolves in an interaction with the 
atmosphere. There is, however, a practical need to start the simulation from an observed fire perimeter, and to modify 
the fire perimeter in a running simulation, which presents a particular problem in a coupled model. The atmospheric 
circulation due to the fire takes time to develop and the heat release from the fire model needs to be gradual, or the 
model will crash due to excessive vertical wind component. 

Therefore, we have proposed creating and replaying an approximate fire history, leading to the desired fire 
perimeter |30|. Replying the fire history allows for graduate release of the combustion heat and allows the atmospheric 
circulation patterns due to the fire to develop. The fire history is encoded as an array of ignition times Ti(x,y), 
prescribed at all fire mesh nodes. To replay the fire, the numerical scheme for advancing O is suspended, and instead 
the level set function is set to 

^(t,x,y) = n(x,y)-t. (5) 

The fuel decay ^ is then computed from Ti, and the resulting heat fluxes are inserted into the atmosphere. After the 
end of the replay period is reached, the numerical scheme of the level set method takes over. 

5. Creating a level set function from a given fire perimeter 

A fire perimeter is considered as a closed curve F, composed of linear segments, and given as a sequence of the 
coordinates of the endpoints of the segments. In practice, such geospatial data are often provided as a GIS shapefile 
ED, or encoded in a KML file, e.g., from LANDFIRE. 

In fSOl, we have proposed a simple scheme for creating an artificial fire history to be used in the fire replay 
scheme ([5]): given an ignition point and ignition time at that point, approximate ignition times Ti at the mesh points 
are established by linear interpolation between the ignition point and the perimeter (Fig.jSj. This simple method was 
already shown to be successful in starting the model from the given perimeter in a simple idealized case (Fig. [4]), with 
the error in the wind speed of only few %. Extensions of the artificial history scheme will be needed for domains 
which are not star- shaped with respect to the ignition point. Running the fire propagation backwards in time to find 
an ignition point is also a possibility, with an intriguing forensic potential |30|. 

The ignition times Ti at locations outside of the given fire perimeter are perhaps best thought of as what the ignition 
times at those locations might be in future as the fire keeps burning. 

Constructing a level set function from a perimeter is one of the basic tasks in level set methods. Given a closed 
curve F, one wishes to construct a function L = L(x,y), such that 

L > outside of F, L < inside of F, L = on F. (6) 
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(a) The difference in the horizontal wind vector (m/s). 



(b) The relative difference in the wind speed. 



Figure 4: The difference in the horizontal wind field at 6.1m between a simulation propagated naturally from a point 
and another one advanced artificially. The first simulation was ignited from a point in the northeast corner of the 
domain, and the fire perimeter was recorded after 40 minutes. This perimeter and ignition location were used to 
generate an artificial history for the first 40 minutes, which was replayed in the second simulation. Both simulations 
were then allowed to advance another 28 minutes. Reproduced from BOl . 



In the application to perimeter ignition, one can then set at a fixed instant t, 

Ti(x,y) = cL(x,y) + it-T), ^{Ux^y) = L(x,y) 

where c is a scaling factor, and proceed with the replay as described in Section |4] 

One commonly used level set function is the signed distance from the given closed curve F, 

L(x,3;) = ±dist((x,3;),F), (7) 

where the sign is taken to be negative inside the region limited by F and positive outside 1321 , and dist stands for 
the Euclidean distance. Surprisingly, such function cannot be defined consistently once the problem is discretized. 
Consider a level set function L that is given by its values on the corners of grid cells, interpolated linearly along 
the grid lines, and F given by its intersection with the grid lines (Fig. [5]). Then, the ratio of the values of L at two 
neighboring mesh corners on the opposite sides of F is fixed by the requirement that L is linear between the two 
corners. In particular, it is not possible in general to define L as the signed distance (|7]). For example, in Fig.[5j the 
ratio L (X) /L (7) is fixed and L (X) does not depend on Z, while L{Y) does. 

One possibility is simply define the values of L next to F by the signed distance, and forget about the exact 
representation of F as L = 0. Instead, in |33 1, we have proposed to find the values of L next to F by least squares. 
Denoting by u the vector of the values L next to F, it is easy to see that u satisfies a homogeneous system of linear 
equations of the form Bu = with at most two nonzeros per rows, and each row corresponding to an edge on the mesh 
that is intersected by F, as the edge XY. We can then find a suitable u minimizing \\u - d\f subject to Bu = 0, where 
d are the signed distances ([7]). Once the values of L near F are found, one can extend L to the whole domain as the 
distance function by the Fast Marching Method (FMM) f34l, or by a simpler and less accurate approximate method 
suggested in |[33]| . 

A better method can be obtained by taking the spread rate into account. The level set function L is a solution of 
the Hamilton- Jacobi equation 



7?|VL| = 1, L = OonF. 
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Figure 5: A level set function linear on the line segments connecting the nodes of the fire mesh cannot be defined at 
the nodes X and Y consistently as the signed distance (|7| from the interface F. The distance of the point X from F 
does not depend on the location of the point Z, while the distance of Y does; yet the values of the level set function at 
X and Y are linear along the segment XY and so fixed by the ratio of their distances from W. 



which can be found by solving the reinitialization equation |[32j Eq. (7.4)] 



dL 



= ±(1-/?|VL|) 



(8) 



where the sign is taken positive outside of F and negative inside. Equation ([8| is solved by upwinding formulas 
moving away from F and starting from the values of L on the other side of F. Alternating the solution process between 
the outside and the inside of F, the values of L on the two sides of F "balance out and a steady- state signed distance 
function is obtained" |[32| p. 66]. 

The situation here is more complicated, because the spread rate R depends on the level set function L following 
([3]). Hence, we freeze L inside R and use successive approximations of the form 



dt 



±\l-R\vi' 



VLk_ 
|VL,| 



|VL. 



'k+i 



6. Data assimilation for the level set fire spread model 

Creating a level set function as in ISOll and in Sec. [5] allows for starting the coupled model from a given fire 
perimeter instead of an ignition point. However, a more general approach is needed for data assimilation. Suppose 
the fire perimeter in the simulation is Fi at time ti. Then at time t2 > ti, the fire evolves to fire perimeter F2. However, 
data is assimilated, changing the state of the fire model and resulting in a diff'erent fire model state with perimeter Fa. 
First we construct level set functions Li, L2, and La for the perimeters Fi, F2, and Fa, respectively, satisfying ([6]). We 
assume that all three level set functions are created using the same method. The resulting approximate formulas will 
be exact in the case of ID propagation with the level set functions linear and having the same slope. They will be 
used point- wise as an approximation otherwise. To emphasize the point-wise application, we write out the arguments 
(x, y) when present. 

6. 1 . Modifying the fire perimeter dynamically 

The state of the atmosphere will no longer match the state of the fire model with the perimeter Fa, and we need to 
make up the evolution of the atmosphere as the fire progresses from the perimeter Fi to the perimeter Fa. Since Fi is 
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completely contained inside Fa, in the region between Fi and Fa, we have Li > and La < 0. The function 

fi,a U, y) = -r-. ^ — —. 7 (9) 

then satisfies 

fi,a = on Fi, < fi^a < 1 between Fi and F^, fi^a = 1 on F2. 

We can then use the function fi^a to create artificial ignition times by 

T, (X, j) = ^1 + ^ , ^'^""^y^ it2 - h) 

Li{x,y)-L^{x,y) 

which interpolates between the perimeters Fi and Fa, and replay the fire history to release the heat into the atmosphere 
gradually, as in Sec.|4] 

6.2. Dynamic estimate of fire spread rate 

A common source of errors in fire modeling is incorrect spread rate. The level set function construction here can 
be used to adjust the spread rate as well. Define /i,2 similarly to ([9]), 

/i,2 (x,y) = —- r — r. (10) 

Li(x,y)-L2(x,y) 

We now use a simple argument of proportions. Assume for the moment ID fire propagation in one direction and that 
fi^a and /i,2 are linear. Then Fi, F2, and Fa are points on the real line, and the spread rates of the simulated fire and 
the spread rate after the data assimilation are, respectively, 

h -ti t2- ti 

However, since /i,2 and /i,^ are linear. 



which gives 



h,2{x)= , 


fl,a (X) = 


x-Ti 

r.-n 


RAx) r^-ri 


/l,2 (X) 


r=tr,. 



R(x) F2-F1 /i,,(x)' " - — ^^^^ 

Thus, using ^ and ([TO]), ^TT) suggests to modify the given spread rate 7? at a point (x, 3;) to become the spread rate 
Ra after the data assimilation, given by 

Li(x,y)-La(x,y) 
RAx.y) = -z—i — r — z— — -R{x,y)- (12) 

Li{x,y)-L2{x,y) 

7. Moisture model 

Fire spread rate depends strongly on the moisture contents of the fuel. In fact, the spread rate drops to zero 
when the moisture reaches the so-called extinction value. For this reason, we have coupled the fire spread model 
with a simple fuel moisture model integrated in SFIRE and run independently at every point of the mesh. The 
temperature and the relative humidity of the air (from the WRF atmosphere model) determine the fuel equilibrium 
moisture contents E 1 35 1, and the actual moisture contents m = m(t) is then modeled by the standard time-lag equation 

dm E - m 

where T^ is the drying lag time. We use the standard model with the fuel consisting of components with 1, 10, and 
100 hour lag time, with the proportions given by the fuel category |36|, and the moisture is tracked in each component 
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separately. During rain, the equilibrium moisture E is replaced by the saturation moisture contents S , and the equation 
is modified to achieve the rain-wetting lag time T^ only asymptotically for heavy rain, 

where r is the rain intensity, tq is the threshold rain intensity below which no perceptible wetting occurs, and r^ is the 
saturation rain intensity, at which 1 - 1/^ ^ 63% of the maximal rain- wetting rate is achieved. The coefficients can 
be calibrated to achieve a similar behavior as accepted empirical models |[37l|38l. See (391I401 for other, much more 
sophisticated models. If moisture measurements are available, they can be ingested in the model ([13] [14]) by a fast and 
cheap Kalman filter in one variable, run at each point independently. 

8. Conclusion 

We have presented new techniques to assimilate the perimeter data at two diff'erent times into a coupled 
atmosphere-fire model, a new method to estimate the adjustment of the model spread rate between the perimeters 
towards the data, and a new coupling of the atmosphere-fire model with a third model, a simple time-lag model of 
fuel moisture. Implementation of the data assimilation is in progress. The moisture model is currently included in the 
code download and will be treated in more detail elsewhere. 
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